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ABSTRACT 

Millisecond pulsars (MSPs) have been studied in detail since their discovery in 1982. 
The integrated pulse profiles of MSPs appear to be stable, which enables precision 
monitoring of the pulse times of arrival (TOAs). However, for individual pulses the 
shape and arrival phase can vary dramatically, which is known as pulse jitter. In 
this paper, we investigate the stability of integrated pulse profiles for 5 MSPs, and 
estimate the amount of jitter for PSR J0437— 4715. We do not detect intrinsic profile 
shape variation based on integration times from ^ 10 to 100 s with the provided 
instrumental sensitivity. For PSR J0437— 4715 we calculate the jitter parameter to 
be fj = 0.067 ± 0.002, and demonstrate that the result is not significantly affected 
by instrumental TOA uncertainties. Jitter noise is also found to be independent of 
observing frequency and bandwidth around 1.4 GHz on frequency scales of < 100 MHz, 
which supports the idea that pulses within narrow frequency scale are equally jittered. 
In addition, we point out that pulse jitter would limit TOA calculation for the timing 
observations with future telescopes like the Square Kilometre Array and the Five 
hundred metre Aperture Spherical Telescope. A quantitative understanding of pulse 
profile stability and the contribution of jitter would enable improved TOA calculations, 
which are essential for the ongoing endeavours in pulsar timing, such as the detection 
of the stochastic gravitational wave background. 
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1 INTRODUCTION 

Millisecond pulsars (MSPs) have be en shown to exhibi t 
highly stable rotational behaviour (e.g. lVerbiest et al]|2009t ). 
They are a vital tool in investigating the physical envi- 
ronment of the pulsars, via precise monitoring of pulse 
times of arrival (TOAs), the so-called pulsar timing tech- 
nique. Previous work using precision pulsar timing has al- 
ready yielded tes ts of General Relativity in the strong 
field regime (e.g. iTavlor fc Weisberd 1 19891 : iKramer et al.l 
I2OO6I ). i mprovements of neutron star equation of state mod- 
els (e.g. iLattimer fc PrakashI l2007l : iDemorest et al.l I2OI0I : 
lOzel et all I2OI0I ). and studies of the interstellar medium 
(ISM) structure (e.g. lYou et "al]|2007l ). Via pulsar timing ar- 
rays (PTAs), upper bounds have already bee n placed on 
the s t ochastic gravitational wave background ( Jenet et al.l 
I2OO6I : Ivan Haasteren et al.l I2OIII : lYardlev et al.|]201ll '). The 
next generation of radio telescopes (such as the Five hundred 
metre Aperture Spherical Telescope, FAST and the Square 
Kilometre Array, SKA) will significantly improve upon the 



available instrumental sensitivity allowing timing to much 
higher precision and for many more sources. This tremen- 
dous advance in hardware will significantly increase the sen- 
sitivity of future PTAs and allow more detailed studies of 
the gravitational wave bac kground, such as its polarisatio n 
and speed of propagation (|Lee et al.ll2008l : iLee et al.ll2010l ). 
Individual gravitational wave sources, such as supermas- 
sive black hole binaries at the centre of g alaxies can also 
be identified via a future PTA study (e.g. iLee et aLll201ll : 
ISesana et aIll2011^ . 

The extreme timing precision required to reach the 
aforementioned scientific goals is more readily achievable 
with MSPs because of their short spin periods, their regular 
rotational behaviour and, last but not least, their greatly 
stable average pulse shapes. In general, single pulses from a 
pulsar show significant shape modulation. Up to now, several 
types of variable behaviour have been observed within the 
population of pulsars. These include intrinsic pulse-to-pulse 
changes caused by random phase jitte r of individual pulses 

systematic position 



cnanges causeg py random pnase ptte r 
(|Cordes fc Downs! 1 19851 : ICordesI [l993h . 
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cha nges of sub-components named "sub-pulse driftins' 



(e.B. Drake fc Craftll 19681 : ICordesll 19751 : lEdwards fc Stapperd 



I2OO3 ). switching between two or more profile shapes on both 
short and long timescales known as "mode-changing" (e. ) 
ICordes et all 1 19781 : iFerguson et"al]|l98ll : iLvne et alJIioiC 
and state cha nges of intrinsic flux density m entioned as 
"nulling" (e.g. iBackeJ Il970l : IWang et al.l 120071 '). Studies of 
bright individual pulses from MSPs have shown only the 
first type of puls e-to-pulse variation (jCognard et al.|[l996l : 
IJenet et al.lll998l ). which can cause fluctuation in the phase 
of integrated profiles and introdu ce TOA uncertainty in ad - 
dition to the radiometer noise ICordes fc Shannon! [201oh . 
This is the effect mainly addressed within the framework of 
this paper. 

Statistically, as the single pulses are clearly unstable, 
shape differences are expected to exist between an integrated 
profile over a short integration time and a standard tem- 
plate formed by averaging many more pulses. Note that the 
shape difference, if sufficiently large, would infiuence the ac- 
curacy of TOA estimation by the standard cross-correlation 
approach (|Liu et al.ll2011^ . Consequently, it is important to 
evaluate whether or not the shape mismatch is significant 
compared with the system noise level. There have already 
been studies investigating the shape correlation between in- 
tegrated profiles (or even single pulses) and a pre-formed 



template (iHelfand et al. 1975 : Rankin fc Rathnasreelll995l : 
iJenet et al.lll998l : iJenet et aLlboOlf ). The results show clear 
shape difference for young pulsars and a less significant or 
even undetectable difference for MSPs with the given sensi- 
tivity. 

If the TOA uncertainties can be shown to be mostly 
due to radiometer noise and phase jitter, based on an as- 
sumed model, the amount of jitter can be estimated from 
timing on a short timescale and such a nalysis has already 
been carried out for PSR J1713-I-0747 (ICordes fc ShannonI 
I2OIOD . The result can be used both for the error analysis 
in timing and as an input for ji t tered shape correctio n ap- 
proaches l|Messenger et al.|[201ll : lOslowski et al.ll201ll '). 

The structure of this paper is as follows. In Section [2] 
we introduce the approach used for evaluating profile stabil- 
ity and jitter estimation. The observations, instruments and 
data pre-processing techniques used are described in Sec- 
tion [3l In Section [4] we present the results before drawing 
our conclusions in Section [S] 



2 PROFILE SHAPE AND PULSE JITTER 
ANALYSIS 

2.1 Stability Analysis 

The shape and phase instability of single pulses will induce 
shape modulation for integrated profiles, which can be mit- 
igated by increasing the number of pulses added. The sim- 
ilarity between an observed integrated profile p and a nor- 
malised standard shape s, obtained from previous observa- 
tions, is simply evaluated by the correlation coefficient p, 
which is defined as: 



(1) 



where i stands for the sample number of the data points 
across the profile. It is shown in Appendix [XI that, assum- 
ing an identical intrinsic shape for the observed profile and 
standard, there is a scaling law between p and the pro- 
file peak signal-to-noise ratio (SNR, defined as pulse peak 
amplitude divided by root-mean-square of the noise) of: 
(1 — p) oc SNR~^. Note that this power law is followed only 
in the high SNR regime (e.g. for SNR> 20). From the scal- 
ing law, we define a shape constant, £, related only to the 
intrinsic profile shape as: 



£ = SNR • ^l-p : 



1/2 



(2) 



where risamp is the number of time samples of a profile and 
i = 1, 2, • • • nsamp. Once the SNR and p are measured, re- 
spectively, the shape constant can then be determined and 
used to compare with the value calculated directly from the 
waveform of the standard. 

If the profile integration increases the SNR as expected 
as SNR oc \/N, where A'' is number of pulses folded, then 
we reach the relation between A'' and p as: (1 — p) oc A'^"^. 
However, the scaling is not obeyed if the SNR varies signifi- 
cantly from pulse to pulse, which can be caused by intrinsic 
fiux variation, system temperature changes and diffractive 
scintillation. So in our analysis we use the effective pulse 
number A'efc, the number of pulses weighted by SNR to ac- 
count for the variation of the profile SNR which d oes indeed 
/Nefc and hence 1 — p oc N^^^ (|Liu et al.l 



scales as SNR oc 
l20Tlh . 



2.2 Phase Jitter 

Single pulse instability can also cause TOA fiuctuations of 
integrated profiles. This phase jitter could, in principle, be 
investigated directly from single pulse data, although the 
sensitivity achieved by the current instruments may not en- 
able sufficient SNR for carrying out the study on all pulses 
within a narrow band. Another approach is to perform tim- 
ing using integrated profiles on short timescales and to es- 
timate the amount of jitter from the timing residuals. As 
a first-order approximation, by assuming an identical shape 
for single pulses and a Gaussian-distributed central phase 
probability, the contribution of pulse j itter to the uncer- 
tainty of TOA can be calculated as in ICordes fc ShannonI 
(|201Cf ). In brief, the measurement error of TO As on short 
timescales (e.g. several hours) can be summarised as: 



2 

""total 



'-' scint 



(3) 



where am, o".t, osdnt and ctq correspond to uncertainty in- 
duced by radiometer noise, pulse jitter, instability of short- 
term diffractive scintillation, and all other possible contribu- 
tions (faults in timing model, instrumental digitisation arte- 
facts, polarisation calibrat i on er r or, etc), respe c tively . Fol- 
lowing iDowns fc ReichlevI (|l983l ). ICordes et all (|l990l ) and 
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ICordes fc ShannonI (|2010l) . we have 



iV X SNRlJ[U'it)Ydt' 



2 



2 

'^scint 



2 / U{t)t^dt 

N JU{t)dt' 



(4) 
(5) 
(6) 



Here A'^ is the number of pulses, SNRi is the signal-to-noise 
ratio for a single pulse, U (t) is the profile waveform, A is the 
sampling time, fj is the width of the Gaussian probability 
of single pulse phase in units of the pulse width, td is the 
pulse broadening timescale, and A'scint is the number of scin- 
tles contained by the integration. We can see that the ratio 
of (Tin to CTj is proportional to the equivalent single pulse 
SNR. The value of SNRi for MSPs, based on current ob- 
servations, is mostly far less than unity, leading to the case 
where the white noise term is dominant in TOA uncertainty. 
However, the contribution by phase jitter will become more 
significant when timing is carried out by the next generation 
of radio telescopes, which will have a significant improve- 
ment in sensitivity of 1 — 2 orders of magnitude over cur rent 
systems (|Nanll2006l : [Schilizzi et alll2007l : [lTu et al.ll20lil ). In 
this case, TOA error prediction based solely on radiometer 
noise will be incorrect. 

The statistics of the timing residuals which are obtained 
by subtracting a timing model from the measured TOAs, can 
be evaluated via a reduced value given by 



Xrec 



N- 



(7) 



where A'^ is the number of residuals, n is the number of 
fitted parameters, and Aa;;, ai are the ith residual and its 
corresponding measurement error, respectively. Normally, ai 
accounts for only the uncertainty of radiometer noise (Jm in 
Eq. Q , w hich can be ob tained from the template matching 
technique (|Tavloij[l992l ). Theoretically, if the timing resid- 
uals are dominated by radiometer noise, a timing solution 
with residuals of Xrcc — 1 will be expected. The existence of 
other types of noise would make the ai underestimated and 
deviate Xiac from unity. The contribution from the diffrac- 
tive scintillation can be estimated by Eq. ((Sjl. Here td can 
be obt ained from the NE2001 Galactic Free electron Density 
Model (jCordes fc Lazicll2002h . and the number of scintles is 
assessable via either a dynamic sp ectrum or more detailed 
calculations in lCordes et al.l (|l990l l. If the additional noise is 
mostly from pulse jitter, one can estimate the jitter noise aj 
by adding its contribution into ai to have Xrec = 1- Accord- 
ingly, the value of the jitter parameter fj can be derived 
based on Eq. ((5]), and the probability density distribution 
of fj can be calculated from the standard distribution, 
given by 



2'=/2r(fc/2) 



fc/2-1 -x/2 

X ' e ' 



(8) 



where x is the value and k — N — n—1 is the degrees of 
freedom. 



Table 1. Details of observations used in this paper. The symbols 
to, A^o and ^obs represent the shortest integration time in the 
dataset, the number of time dumps and the observational length, 
respectively. The asterisk (*) denotes that the data were pre- 
processed by the DFB and the others are from CPSR2. 



Pulsar Name 


MJD 


Receiver 


to (s) 




Tobs (hr) 


J0437-4715 


53576 


MB 


16.7772 


1595 


8.7 




53621 


MB 


16.7772 


1498 


8.9 




53864 


MB 


67.1088 


261 


6.9 




54095* 


H-OH 


60.0017 


125 


2.4 




54222 


H-OH 


67.1088 


135 


9.8 




54226 


H-OH 


67.1088 


152 


8.9 


J1022-I-1001 


53260 


MB 


16.7772 


114 


0.5 


J1603-7202 


53166 


H-OH 


16.8099 


114 


0.5 


J1713-I-0747 


53221 


H-OH 


16.7771 


215 


3.3 


J1730-2304 


53145 


H-OH 


16.7772 


114 


0.5 



3 OBSERVATIONS 

Frequent observations, typically weekly, of more than 
20 MSPs are performed at the Parkes 64-m radio 
telescope (|Verbiest et al.l |2009| ). Here we use the data 
from five sources (PSR J0437-4715, PSR J1022-fl001, 
PSR J1603-7202, PSR J1713-f0747, PSR J1730-2304), 
collected using either the Parkes 20 -cm multibeam (MB) 
receiver (|Stavelev-Smith et all Il996l ) or the 'H-OH' re- 
ceiver. The data were processed online with the Caltech- 
Parkes-Swinburne Recorder 2 (CPSR2), a 2-bit coherent de- 
disperser back-end th at records two 64-M Hz wide observing 
bands simultaneously (|Hotan et al.ll200^ . These bands were 
centred at observing frequencies of 1341 and 1405 MHz, re- 
spectively. Data collected before MJDs 53740 were folded 
every 16.7772 s, and after that every 67.1088 s. Off-source ob- 
servations of a pulsed noise probe at 45° to the feed probes, 
but with otherwise identical set-up, were taken at regular 
intervals for the purpose of polarimetric calibration. Addi- 
tionally, we analysed data for PSR J0437— 4715 from the new 
Parkes Digital Filterbank (DFB) system, a digital polyphase 
filterbank capable of 8-bit sampling. In Table [1] we present 
the details for all selected datasets. 

The data were then pre-pro c essed with the psrchive 
software package (| Hot an et al.l |2004| ) . Specifically, for 
CPSR2 data we correct ed the 2-bit digitisation artefact 
(|jenet fc Andersonll 19981 ) by applying the algorithm in van 
Straten (in preparation), and removed 12.5 % of each edge of 
the bandpass to avoid possible effects of aliasing and spectral 
leakage. A full receiver model was determined and applied 
to pe rform the polarisa tion calibration to the MB receiver 
data (Ivan Straten|[2004[ ). as the receiver suffers from strong 
cross-coupling between the orthogonal feeds. For the H-OH 
data from both back-ends we used the common single-axis 
model instead, as the co upling was found to be an order 
of magnitude weaker (e.g. [Manchester et al.|[2"010l ). The sig- 
nals from each polarisation were summed into total power 
(Stokes /), while 0.5 MHz frequency channels were kept for 
later analysis purposes. The template profiles used for the 
correlation coefficient calc ulation and t he cross-correlation 
procedure to estimate CTrn (|Tavloilll992l ). were obtained in- 
dependently from the datasets shown in Table [1] All CPSR2 
datasets have passed through the test shown in I Liu et al] 
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l|201ll ) to ensure no significant 2-bit digitisation distortion 
is present. The test also showed that the template matching 
produced the radiometer noise uncertainty as expected by 
Eq. @. 



4 RESULTS 

4.1 Measurement of shape constant 

The measurement of the shape constant for individual in- 
tegrations gives an estimate of the profile stability, which 
here we have carried out for the aforementioned five MSPs. 
For PSR J0437-4715, £ was determined using the the 
first 500 time dumps of the MJD 53621 dataset. For 
PSRs J1022-I-1001 and J1730-2304 we integrated the time 
dumps into 1.0 and 1.8 minutes, separately, so as to obtain 
sufficiently high SNR for calculation of £ and statistical anal- 
ysis (see e.g. Fig. lAlf) . The p is computed with respect to 
the on-pulse phase and the root-mean-square (RMS) devia- 
tion of the noise is estimated based on the off-pulse region, 
from which one can derive the value of £ from Eq. The 
expected values of £ (denoted by £o) are calculated directly 
from the shape of the high SNR standards formed from pre- 
vious observations, which are also used in the calculations of 
p. The errors in p and the RMS of the baseline of the profile 
are given by Eq. (|A6|) and by ~ RMS/y27isamp, respectively. 

We present the shape constant measurement of 
PSR J0437-4715 for the MJD 53621 dataset in Fig. [T] as 
an example, and show the statistical result in Table [2] It is 
clear that for most sources, within the range of estimated 
error the measured £ is in accordance with the analytical 
value. The RMS of the measured value also matches the 
mean error bar for each data point. The deviation of the 
measured £ from the expectation for PSR 1730—2304 may 
indicate an intrinsic level of profile shape change, but could 
also be due to an insufficient number of measurements. Note 
that, in this case, the statistics are only based on 13 data 
points, while for each of the others more than 30 measure- 
ments were available. 

It is worth noting that there has already been pre- 
vious work regarding the stability of PSR J1022-I-1001 
on bo th short and long t imescales, with inconclusive re- 
sults ([Kramer et al. I ll999l : iRamachandran fc Kraiiieil l2003l : 
iHotan et al.ll2004 ). On timescales of roughly an hour, our 
results do not show profile changes. However, the analysis 
was carried out using the whole on-pulse phase, and may not 
be sensitive to variations occurring only around the profile 
peaks. Any shape variations, if they exist, are likely to affect 
mostly the peaks. A detailed analysis of the statistics of the 
peak ratios will be published by Purver et al. (in prepara- 
tion). 

To investigate the scaling relation indicated by Eq. ([2]), 
for each source we perform incremental integrations along 
the datciset, adding more pulse periods each time. Our re- 
sults are presented in Fig. [2] which shows the relation be- 
tween the weighted number of integrated pulses and p. It 
is clear that all curves are linear in log-log space for the 
Nefc — (1 — p) relation as expected. The fitted slopes all lie 
in the range of (0.94, 1.00), as given in the last column of 
Table [5] This result, together with the £ statistics in Ta- 
ble O indicates no detectable profile shape variation along 
the integration. 




-0.04 -0.03 -0.02 -0.01 0.01 0.02 0.03 0.04 

MJD-53620.7 

Figure 1. Shape constant measurements for the first 500 inte- 
grations of the MJD 53621 observation for PSR J0437-4715. The 
time baseline of these measurements is ~ 1.9 hours. 



Table 2. Statistical results of the profile shape constant mea- 
surement for all five sources. The parameters tj^t, £o> RMS, 
a and Ainc represent integration time of the integrated profiles, 
shape constant calculated from the standard, mean of measured 
C, RMS of measurements, mean of individual measurement er- 
rors, and the fitted slopes of N^f^ versus (1 — p) curves in Fig.O 
respectively. 



Pulsar Name tint (s) Co 



RMS 



J0437-4715 


16.8 


4.8 


4.9 


0.32 


0.31 


-1.00(1) 


J1022-I-1001 


67.1 


2.7 


2.6 


0.22 


0.20 


-0.98(1) 


J1603-7202 


16.8 


3.9 


4.0 


0.33 


0.30 


-0.99(1) 


J1713-I-0747 


16.8 


3.7 


3.8 


0.23 


0.22 


-0.94(2) 


J1730-2304 


117 


2.9 


3.2 


0.24 


0.27 


-0.98(1) 



4.2 Fit of phase jitter 

Estimation of the jitter parameter fj can be performed 
only on the brightest source PSR J0437— 4715, as single 
pulse SNR for the other sources is not high enough and 
radiometer noise is still dominant in the timing residuals. 
Here we use all PSR J0437-4715 datasets listed in Table [H 
each of which contains over a hundred individual integra- 
tions to yield a stable statistical result. When performing 
the timi ng of the datas e ts, w e used the timing models de- 
rived bv IVerbiest et al.l (|2009l ). TOAs and their uncertain- 
ties were determined through cross-correlation with a pre- 
formed standard, individually for each side of band. No 
EFAC or EQUAD value were applied to change the measure- 
ment precision as done c ommonly in timing analysis (e.g. 
Ivan Haasteren et al.ll201lf l. 

The timing solution of PSR J0437— 4715 shows a widely 
scattered series of TOAs over a timescale of several hours, 
with a reduced fa-r larger than unity, which strongly in- 
dicates the existence of extra uncertainty contributions be- 
sides radiometer noise. Fig. [3] shows the timing residuals of 
the MJD 53621 dataset for the 1405 MHz band. The TOAs 
were produced by summing over the whole bandwidth of ef- 
fectively 48 MHz, yielding reduced values of 8.45 and 8.46 
for the 1341 MHz and the 1405 MHz side band, respectively. 
Note that from Eq. ([6]), the TO A fiuctuation induced by un- 
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1000 10000 100000 le+06 



Nefc 

Figure 2. N^ic — (1 — p) plot for all five sources. The effective 
pulse number is used to compensate for the varied weight of each 
integration due to the difference in SNR. All curves show expected 
linear behavior with fitted slopes to —1 in log- log space. 




-0,2 -0,1 0,1 0,2 

MJD-53620,9 



Figure 3. The short-term timing solution of PSR J0437-4715 
over ~8 hours dataset on MJD 53621 at a central frequency of 
1405 MHz with 48 MHz bandwidth. 

stable profile broadening due to the ISM is approximately 
the scattering timescale once the observational bandwidth is 
much narrower than the scintillation frequency scale (thus 
A^'acint ~ 1). For PSR J0437-4715, as the scintillation band- 
width and pulse broad ening time are ~ 0.4 GHz and ~ 1 ns 
l|Cordes fc Lazioll2002l '). respectively, this uncertainty is neg- 
ligible for the datasets used in this paper. 

Following the method mentioned in Section 12.21 we 
measure fj based on the shortest integrations from all 
PSR J0437-4715 datasets in Table [T] and the results are 
summarised in Table [S] It is clear that the measurements 
from datasets collected by different receivers are consistent 
with each other. This does not indicate significant uncer- 
tainty contribution by polarimetric calibration error. Re- 
sults from two types of back-end also show consistency af- 
ter correcting for the bias induced by low-frequency noise 
of the CPSR2 data (see the next paragraph for details). 
Additionally, the same result for fits from both sides of 
the band suggests that the intensity of phase jitter does 
not vary significantly on small frequency scales (there is a 
64 MHz difference between the two bands) at observing fre- 



Table 3. Results of jitter parameter fj measurements and K- 
S tests of all PSR J0437-4715 datasets from two 48-MHz sub- 
bands. The parameters / and V represent the central frequency 
and the p-value of K-S test, respectively. 



Dataset (MJD) 


/ (MHz) 




fj 




V 


53576 


1341 





067 ±0 


004 





89 




1450 





067 ±0 


004 





93 


53621 


1341 





067 ±0 


004 





98 




1405 





066 ±0 


004 





62 


53964 


1341 





066 ±0 


004 





88 




1405 





065 ±0 


004 





92 


54095 


1341 





069 ±0 


006 





90 




1405 





073 ±0 


006 





92 


54222 


1341 





069 ±0 


007 





52 




1405 





067 ±0 


007 





77 


54226 


1341 





068 ±0 


005 





20 




1405 





066 ±0 


006 





64 



quencies around 1.4 GHz. The combination of all measure- 
ments for both sides of the band achieves an estimated fj 
of 0.067 ± 0.002. To study the statistics of the residuals, we 
perform Kolmogorov-Smirnov (K-S) tests on each dataset, 
on TOAs weighted by the modified uncertaint ies account- 
ing fo r phase jitter. The measured p- values (e.g. lPress et al.l 
Il986t ) are summarised in Table [3l not suggesting a signif- 
icant deviation of the weighted residuals from a Gaussian 
distribution. This demonstrates the dominance of Gaussian 
noise in the residuals. 

To investigate the dependence of the result on the 
length of individual integrations, we calculate the weighted 
RMS against integration time tint for multiple combina- 
tions of instruments. In detail, MJD 53576 (CPSR2-fMB), 
MJD 54095 (DFB-hH-OH) and MJD 54226 (CPSR2+H- 
OH) datasets are chosen to demonstrate different hardware 
configurations, and the results are shown in Fig. U Clearly, 
the curves yielded by CPSR2 data collected by two receivers 
are consistent, again implying no significant contribution of 
timing uncertainty by polarisation calibration. The relation 
obtained from DFB data achieves a fitted slope close to 
—0.5, supporting the idea from Eq. ([3} that timing resid- 
ual scales as the square-root of the number of integrated 
pulses once only the uncertainties by radiometer noise and 
jitter are significant. The residuals from the CPSR2 data co- 
incide with those from the DFB in 1-min integrations, and 
then saturate at the level of «150ns as the integration time 
is extended. This implies that the difference may be due 
to different digitisation procedures. The saturation corre- 
sponds to additional self-correlated noise which contributes 
«8% bias in our fj measurement based on 1-min integra- 
tions of CPSR2 data, which has been corrected for in the 
results of Table |3l 

The expression for ctj indicates that the uncertainty due 
to jitter is independent of the observational bandwidth and 
central frequency. To illustrate this, the MJD 53621 dataset 
is divided into two, three, four and six sub-bands each time. 
Then a fit for the jitter parameter is carried out individually 
on each sub-band and the results are combined incoherently 
to obtain an estimated value of fj for a given bandwidth. 
The procedure is carried out on both side-bands and the 
result is shown in Fig. [5] It is clear that at both central fre- 
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Figure 4. Weighted timing residual versus integration time re- 
lation for three combinations of instruments. The curve yielded 
by the DFB data shows the expected power law with an index 
close to —0.5, while the curves from CPSR2 data clearly deviate 
from this, which implies additional uncertainty perhaps duo to 
the low-bit digitisation procedure. 
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Figure 5. Fitted amount of jitter against chosen bandwidth for 
both sides of band from MJD 53621 dataset. No clear dependency 
between these two parameters is visible. 
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Figure 6. Estimated jitter parameter based on an 8-MHz sub- 
bands from the MJD 53621 dataset. The fj shows no clear trend 
of evolution across a frequency range of ~ 100 MHz. 



5 CONCLUSIONS AND DISCUSSIONS 
5.1 Summary of the results 

In this paper, we investigate the issue of MSP profile sta- 
bility based on five pulsars in total. A shape constant asso- 
ciated with the correlation coefficient is defined to quantify 
the stability. No significant shape modulation of integrated 
profiles beyond the measurement error is found for integra- 
tion times from 10^ to 10^ s. For PSR J0437-4715 we es- 
timate the jitter parameter by performing timing on short 
timescales and comparing the actual timing residual with 
the amount expected from radiometer noise. The fitted fj 
is found to be consistent on both sides of the bands (64 MHz 
difference in central frequency) , and the combination of sev- 
eral datasets results in an estimate of 0.067 ± 0.002. It is 
also demonstrated that all the other potential sources of 
TOA uncertainty, besides radiometer and jitter noise, do not 
strongly influence the measurement. Additionally, we show 
that jitter noise scales neither with bandwidth within a 50- 
MHz band nor with frequency across a range of ~ 100 MHz 
at 1.4 GHz, which supports the idea that pulses are equally 
jittered on small frequency scales. Such results, if still valid 
for wider frequency range, would suggest that the jittering 
uncertainty cannot be mitigated by exten ding the observing 
bandwidth (also see lOsIowski et al]|201ll ). 



quencies jitter noise remains once the bandwidth is changed. 
The result indicates that, on small frequency scales, pulses 
are jittered in the same way, otherwise jitter noise would 
not remain the same after summing the entire bandwidth. 
In Fig. |S] we plot the fitted fj based on 8 MHz bandwidth 
against the central frequency for each sub-band. The group 
of fj values do not show a clear dependence on the observa- 
tional frequency, and yield average, weighted RMS and re- 
duced of 0.069, 0.006 and 1.1, respectively. The mean cor- 
relation coefficient between residuals of different sub-bands 
were calculated to be ~ 0.51, which indicates a correlation 
among the time series and supports the idea of equal jitter- 
ing on small frequency scales. 



5.2 Future telescopes 

Apart from bright single pulses and giant pulses (e.g. 
ICognard eta\] Il996l : IJenet et al1ll998l '). pulse jitter of the 
majority of MSP pulses is still not detectable with currently 
available sensitivity. However, with the next generation of 
radio telescopes, the shape modulation of integrated pro- 
files for some of the bright MSPs will become visible. In 
Fig. [T] based on the aforementioned jitter model and assum- 
ing fj — 0.1, we perform a Monte Carlo simulation to cal- 
culate the relation between the number of integrated pulses 
and 1 — p for the case of jitter only, and compare the re- 
sult with the curves from considering radiometer noise only 
for a few instruments. We assume a 5-ms period, a 100- 
fis pulse width, and a 5-mJy flux density at 1.4 GHz for a 
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Figure 7. 1 — p versus A'^ curve prediction of PSR J0437— 4715 
for different instruments and simulation. The pulses are created 
with identical shape and a jitter parameter of fj = 0.1 is used to 
simulate the Gaussian-distributed phase variation. 

sample MSP, and 1.4 GHz frequency with 300 MHz band- 
width for observation. The gains of FAST and S KA are as- 
sumed to be 20 m^/K and 100 m^/K, respectively (jNanl2006l : 
ISchilizzi et alll2007l '). For the assumed jitter model the scal- 
ing also follows 1 — p oc A'^"^ as calculated in Appendix iBl 
It is shown that an SKA observation of an MSP of typical 
brightness, pulse jitter is comparable to radiometer noise 
in influencing the correlation-coefficient value. Note that in 
the applied jitter model all single pulses are assumed to be 
identical, so the simulated jitter curve can potentially move 
upwards once the shape modulation is also accounted for. 
Future observations with the SKA of PSR J0437-4715 will 
be totally dominated by pulse jitter in shape variation, so it 
will be ideal for single pulse study. 

Once the shape variation of integrated profiles by 
pulse jitter becomes significant, the current cross-correlation 
method for the measurement of TOAs would fail in estimat- 
ing the TOA uncertainty correctly. Specifically, the model 
in the template-m atching procedure now needs to be of the 
form (|Tavloilll992l ): 

Pit) = Ao*S{t) + n{t)+j{t), (9) 

where P{t) is the observed profile, S{t) is the template, and 
n{t) is the noise function. The parameter j{t) is the jitter- 
induced shape perturbation, which, referring to MSPs, is 
mostly negligible compared with n{t) for the current sensi- 
tivity. If this shape difference becomes sufficiently large, the 
calculated T OA precision wo uld fail to follow the expected 
SNR scaling (jLiu et al.ll2oTl] ). In this case, the shape modu- 
lation would need to be modelled (e.g., by principal compo- 
nent analysis, see ICordes fc ShannonI I2OI0I : lOslowski et al.l 
I2OIII ) and then a global determination of the unknown pa- 
ram eters could be performe d to properly estimate the TOA 
(e.g. [Messenger et al.ll201ll '). 
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APPENDIX A: CORRELATION-COEFFICIENT 
SCALING FOR ADDITIVE NOISE 

Assume the observed profile is a superposition of a nor- 
malised template p and Gaussian noises n, i.e. pi = Si +ni. 
The subscript i goes from 1 to n, the number of the profile 
bins. We regard p and n as n-dimensional vectors. Follow- 
ing Eq. ([T]), the correlation coefficient between the perfect 
template s and observed profile p is 



P = 



Ci(ci -I- Hi) 

i 



(Al) 



where d = Si — s. 

Assume that the noise n is a multivariate Gaussian 
with probability distribution /(n) and covariance matrix C, 
where dj = a^Sij. The expectation value of correlation co- 
efficient (p) and its second-order moment (p^) are 



(P) 



pf{n) dn 



P /(n) dn 



(A2) 
(A3) 



O 
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1 2 
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Figure Al. Theoretical and simulated results presenting the 
relation between SNR and correlation factors. In each case the 
dashed line represents the analytic solution and the solid line 
stands for the numerical result. The parameter nbin is the num- 
ber of segments for the profiles. 



The variance for p is then af, = (p ) — (p) . The 
Eq. (|A2|I and HA3P can be i ntegrated by transforming into 
hyper-spherical coordinates (|Mathews fc Walked 1 1970 ). Al- 
though no analytical expression could be found, by using 
asymptotic technique we derive the results for the case of 
a large sample number which match both the low and high 
SNR cases as below: 



l-{p) 



X 



1/2 



X< 1 



1 1/2 



nsamp(l + X) 



X> 1 



(A4) 
(A5) 

(A6) 



where 



(A7) 



In Fig. lAll a set of Monte Carlo simulations are per- 
formed so as to test the validity of the derivation. Here a 
simple Gaussian template shape is assumed and profiles are 
created by adding normal distributed noise onto the tem- 
plate. It can be seen that for most SNR ranges the results 
from both approaches coincide with each other, and for high 
SNR value, 1 — p scales linearly with the increase of signal. 
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APPENDIX B: CORRELATION-COEFFICIENT 
SCALING FOR JITTER 

In this section, we prove that the relation 1 — (p) oc 
Assuming single pulses are of identical waveform po [cj>) and 
the effect of phase jitter is to introduce a random phase 
to each single pulse, i.e. the i-th single pulse takes a wave- 
form of Po [4> + 1 where the A</)i is a random phase. The 
waveform of the template profile p((/>) is defined by summing 
infinite number of single pulses as 



1 



(Bl) 



Meanwhile, an integrated profile pn obtained from an 
observational session is yielded by averaging A'' single pulses 
(iV > 1) as 



1 ^ 



(B2) 




Figure Bl. The 1 — (p) as function of N and crj, where ctj is 
the variance of the phase jitter probability density function. The 
dashed lines are from analytical calculation, the Eq. I IBIOII . while 
the solid lines are from direct Monte Carlo simulations. 



The correlation coefficient between the template p{cj)) and 
A''-averaged integrated pulse profile Pn{4i) is then 



p{<t>)PN{(t>) dct> 



p{(t>f d4> I pN{4>)^d(f) 



(B3) 



The evaluation for the statistical e xpectation of corre- 
lation coefficient is already presented in lCordes fc ShannonI 
l|2010l ). In this appendix, a slightly different approach is ap- 
plied. Letting 



we have 



/o/(0)+P(0M(0)rf<^ 



[/; p2(0) # pH^) + 2p{<i>)5{^) + 52 (^) ■ 



(B4) 



(B5) 



Since the number of pulse in observations is usually large, 
we expect that the PN{<i>) is very close to p((^), which leads 
to S[(t)) <^ p(0). In this case, one can show that 



2\V loP'Wd^ ) loP'Wd^l U 



(B6) 



which is equivalent to 



and 



1 \ 2 

p{4>)PNi<t>) d(j) 



7V2 



p {(f>)d(f> 



+ ^ "^^l' #'p(</')p(0'){Po(<^ + A0)po(0' + A^)) . (B9) 

Here we assume that the phase jitter is independent, i.e. 
(A0»A0j) = 0, if i / j. Thus 



(BIO) 



where 
5 = 

K = 



/;{pg(^+A0))d0 



IoP'i4>)d4> ' 

dcbj^l d,^'p(<^)p(0'){Po(0 + A,^)po(<^' + A,^)) 



(Bll) 



(/,>2(0)d,^)' 



(BI2) 



Because neither S nor K is A^-dependant, the Eq. (|BIO|l 
clearly show the scaling relation that 1 — (/o) oc 1/N. 

We also further compare the result of f — (p) form Monte 
Carlo simulations and from the Eq. (|Bf OP in the Fig. IBll 
Clearly, the numerical simulation and analytical calculation 
match each other at larger A*" limit as we expected. 



, l/ loP%Wd^ ( j^p{^)pN{^)d 

^\I^pH4>)d4> V /oP'('^)rf'^ 

■(B7) 

To understand how the 1 — (p) depends on the number of 
pulses A*', we have t o expand PN{<t>) in the ab ove equations. 
One can then have (jCordes fc Shannonll201(il ) 



PnW d(t> 



- N 



+ 



iV2 



N 



(p^(0 + A0))#, 



(B8) 



